home *** CD-ROM | disk | FTP | other *** search
/ SPACE 1 / SPACE - Library 1 - Volume 1.iso / program / 441 / fplib20 / dnotes < prev    next >
Text File  |  1990-11-23  |  5KB  |  134 lines
  1.         Design notes for FPLIB Release 2.0
  2.         ==================================
  3.  
  4.  
  5. NUMBER REPRESENTATION
  6. ---------------------
  7.  
  8. The floating point format used is Motorola's "Fast Floating Point":
  9.  
  10.         -----------------------------------------
  11.         |    24 bits mantissa    |S| 7 bits exp |
  12.         -----------------------------------------
  13.  
  14. The mantissa has an implicit binary point just off to the left.  The
  15. exponent is binary, stored excess-64.  In other words, the actual
  16. mathematical exponent has 64 added, and the resulting 7-bit number is
  17. unsigned.
  18.  
  19. Negative numbers are stored as sign-and-magnitude.
  20.  
  21. All floating-point operations expect normalized numbers, and generate
  22. normalized numbers if the inputs are normalized.  A number is
  23. normalized if the most significant bit is 1 (except for 0.0, which is
  24. stored as 32 bits of zero, never negative).  Also, a nonzero mantissa
  25. with a biased exponent of 0 is not normalized, letting us test the low
  26. byte or whole long as a test for 0.0.  This all means that the most
  27. significant bit is actually redundant, and could be left off, giving
  28. one extra bit of precision.  Oh well.
  29.  
  30. Some of these routines make use of the assumption that all arguments
  31. are normalized, and may generate slightly bogus results if that isn't
  32. true.  Release 1.0 assumed that a biased exponent of zero could have
  33. a nonzero mantissa; that has been fixed.
  34.  
  35. Some examples: 1.0 is hex 80000041 , or 0.5 * 2 ** 1 (binary 0.1,
  36. exponent 1 + bias of 64).  -0.75 is hex C00000C0 (note the sign bit is
  37. inverted, but no other bit is).
  38.  
  39.  
  40. ACCURACY AND LIMITS
  41. -------------------
  42.  
  43. The accuracy is in the 24-bit mantissa, which gives between 7.22 and
  44. 7.52 decimal digits of relative error.  Decimal constants are given to
  45. 8 digits, and the compiler is trusted to convert them right (it
  46. generally does a reasonable job).
  47.  
  48. The range is in the exponent.  The smallest nonzero and largest
  49. positive normalized numbers are ~5.4e-20 and ~9.2e18.
  50.  
  51. Rounding of ambiguous numbers uses round-to-even (see D. E. Knuth's
  52. Seminumerical Algorithms for a discussion on this).  This means that
  53. if the trailing portion is exactly one half of the least significant
  54. bit, the lsb is forced to be the nearer even number.  We often keep
  55. several insignificant bits to be sure this is done correctly.
  56.  
  57.  
  58. CODING TRICKS
  59. -------------
  60.  
  61. It's quite normal in a library of this type to resort to bit-twiddling
  62. in the internal representation.  To help this, FPLIB.H defines a C
  63. union, thus:
  64.  
  65. union    {
  66.     unsigned long ul;
  67.     float f;
  68.     char sc[4];
  69.     };
  70.  
  71. The number can be looked at as a floating or integer long (for handling
  72. bits in the mantissa), the sign is the sign of sc[3], and the exponent
  73. is the rest of sc[3].  Other shortcuts are commented in the inividual
  74. routines.
  75.  
  76.  
  77. MATHEMATICAL INFORMATION
  78. ------------------------
  79.  
  80. Well, this is mostly what I remember from my undergraduate days.  Many
  81. of these functions have a precise mathematical formula or algorithm
  82. that will give the result, but is too slow for practical use.  Instead,
  83. we look for an approximation that is quick to compute, and gives us
  84. just enough precision and no more.  After all, why compute to 100
  85. digits accuracy when you can only represent less than 8?  Some of these
  86. approximations look very strange, but they all have sound mathematical
  87. bases.  The name Chebyshev springs to mind.
  88.  
  89. Dedicated mathematicians have struggled for decades to bring you these
  90. formulae and, just as important, a predictable accuracy of better than
  91. 7.5 decimal digits.  Because of this, we can execute a known number of
  92. operations instead of loop-and-test (as you would in, for example, a
  93. simple Newton's iteration for square root).
  94.  
  95. The approximations are generally polynomials, usually avoid division, and
  96. are valid only over a small range of the input parameter.  We use simple
  97. mathematical identities to extend the range.  For example, the cyclic
  98. nature of the trigonometric functions means we only need calculate in the
  99. range 0..pi/2.  For another example, we can do much of the square root
  100. calculation by the simple expedient of dividing the binary exponent by 2.
  101.  
  102.  
  103. RANDOM LIBRARY INFORMATION
  104. --------------------------
  105.  
  106. Release 1.0 of Sozobon C had a few problems in the floating support:
  107.  
  108. - Two negative floating numbers compared wrong (returned > when it
  109.   should be <).  I corrected that.  However, the other routines don't
  110.   assume the correction has been installed in the library (which it
  111.   wasn't, at first...).
  112.  
  113. - The library was not optimally ordered, and had to be scanned twice.
  114.   I straightened that out.
  115.  
  116. - Add and subtract were truncating.  Rounding is generally a better idea.
  117.  
  118. atof() calls sscanf().  Ick.
  119.  
  120. -------------------------------------------------------------------------------
  121.  
  122. These routines may be used without restriction, but I retain the copyright
  123. on them, and the Only True Upgrades will be those I make.  If you have any
  124. improvements, please send them to me.
  125.  
  126. March 1989
  127.  
  128. David W. Brooks
  129. 36 Elm Street
  130. Medfield, MA 02052
  131.  
  132. BROOKS@CSSS-A.PRIME.COM
  133. {uunet,mit-eddie}!csss-a.prime.com!brooks
  134.